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Abstract 

We study the temperature dependence of discretization errors in nuclear lattice simulations. We 
find that for systems with strong attractive interactions the predominant error arises from the 
breaking of Galilean invariance. We propose a local "well-tempered" lattice action which eliminates 
much of this error. The well-tempered action can be readily implemented in lattice simulations 
for nuclear systems as well as cold atomic Fermi systems. 
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I. INTRODUCTION 



Nuclear lattice simulations address the nuclear many-body problem by combining numer- 
ical lattice methods with effective field theory. There have been several recent studies on 
the subject of nuclear lattice simulations , [lO| • The starting point 

is the usual starting point of effective field theory. All local interactions consistent with 
the symmetries of low-energy nuclear physics are organized by counting factors of Q/Ahigh, 
where Q is the typical nucleon momentum scale and Ahigh is the high-momentum scale 
where the effective theory eventually breaks down. For chiral effective field theory we have 
Ahigh ~ 47r/^, and for effective field theory without pions Ahigh ~ "^tt- Interactions in the 
effective theory are truncated at some order in Q/Ahigh? and the remaining interactions are 
put on a space-time lattice. Coefficients for the interactions are determined by matching 
to scattering phase shifts and few-nucleon spectra. Once the interaction coefficients of the 
lattice effective theory are determined, the many-body system can be simulated nonpertur- 
batively using Monte Carlo. The method can be applied to nuclei at zero temperature as 
well as the thermodynamics of nuclear and neutron matter at nonzero temperature. Similar 
lattice effective field theory techniques have been used to study cold atomic Fermi systems 
in the limit of short range interactions and long scattering length [ll, 
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In addition to Q and Ahigh there is also a momentum cutoff scale A. On the lattice with 
lattice spacing a, the cutoff momentum scale is A = vra"^. Ideally one should increase A 
systematically to extrapolate to the continuum limit A — oo. For any finite set of diagrams 
with the required local counterterms this is not a problem. However when diagrams are 
iterated to all orders complications arise when the interactions involve singular potentials. 
For example in pionless effective field theory it is known that a three-body counterterm is 



required at leading order 
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24 1251 . 126( 1 . With the three-body countert- 



erm in place the continuum limit is well defined for few-body calculations. Unfortunately at 
very high cutoff momentum this approach involves removing spurious deeply-bound states 
by hand, and there is no way to do this in many-body simulations. The problem is no 
better in effective field theory with pions. In this case the pion tensor force generates in- 
stabilities in higher partial wave channels at large A 
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301 ]. In short the presence 



of continuum limit instabilities and computational constraints means that for many-body 
simulations one is restricted to a finite range of values for A. Therefore it is important to 
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understand and control errors that occur at finite cutoff momentum. 

In tliis study we investigate finite cutoff errors on the lattice at nonzero temperature. 
While our results and conclusions apply to general few- and many-body nuclear systems, 
we center our discussion on simulations of dilute neutron matter. In particular we consider 
an idealized limit of neutron matter with zero range two-body interactions. We take this 
zero-range two-body contact interaction as the only interaction. Therefore Ahigh oo, and 
it is straightforward to show that this idealized theory has no continuum limit instabilities. 
For these reasons it is a useful testing ground to study Q/A cutoff errors without additional 
complications. Zero-range neutron matter is a good approximation to actual dilute neutron 
matter when the spacing between neutrons is sufficiently large. This occurs at about 1% of 
normal nuclear matter density or less. 

Our interest in finite cutoff errors at nonzero temperature is motivated by a recent anal- 
ysis of zero-range neutron matter on the lattice which found sizable lattice errors at nonzero 
temperature Tjl. This stands in contrast with zero temperature simulations which found 

nn 

little dependence on lattice spacing [IJ, 117|. The lattice spacing dependence at nonzero 
temperature was first noticed in the results of many-body lattice simulations and then ana- 
lyzed by calculating coefficients of the virial expansion. The second-order virial coefficient 
b2(T), where T is temperature, was found to be too large when computed on the lattice. 
While the source of the error was unknown, it was suggested that tuning the two-body 
interaction to give the correct value for b2{T) might improve the results of the many-body 
simulation. This suggestion was carried out in [8|, and the many-body lattice results with 
the retuned interaction showed little residual dependence on lattice spacing. Similar cutoff 
errors were found in 15|, ll6|. However the analysis did not distinguish between cutoff errors 
due to nonzero temperature and cutoff errors due to nonzero density. 

In this paper we answer some of the questions raised by the findings in [3] and Q]. In 
particular we discuss the source of the large temperatiire-dependent lattice errors, why the 
measured energies tended to be too low, and why in [8] it was possible to cancel much of 
the error by retuning the two-body interaction. We also propose a simple modified lattice 
action which eliminates most of the large temperature-dependent lattice errors from the 
beginning. The results of our analysis should be useful for reducing systematic errors in 
future nuclear lattice simulations as well as other strongly-attractive fermionic systems. 
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II. VIRIAL EXPANSION 



The virial expansion for the equation of state has been used to st udy neutron and nuclear 



matter as well as fermionic atoms near the classical regime [31|, 123, 133|, |3J] The virial 
expansion can be regarded as a power series in fugacity, z = e^^, where j3 is the inverse 
temperature and /i is the chemical potential. For example the logarithm of the grand 
canonical partition function per unit volume for neutron matter can be written as 

Ur,ZG=^[z + b,{T)z^ + h{T)z' ■ ■ ■] , (1) 



where 



A. - J'-^ (2) 



is the thermal wavelength and m is the neutron mass. We can use the virial expansion 
to compute thermodynamic observables when the thermal wavelength is smaller than the 
interparticle spacing. The neutron density, p, can be computed by taking a derivative of 
In Zg with respect to the chemical potential, 

To second order in the virial expansion we find 

p = ^[z + 2h{T)z' + ...]. (4) 
Taking into account Fermi statistics, we get 

^free^y) = -2'^/'^ ^ -0.177 (5) 

for a free gas of neutrons. 

With the interactions turned on, the second virial coefficient can be computed by ex- 
tracting the term in the partition function proportional to z'^, 

hiT) - b'riT) = ^ {Tr,[exp{-(3H)] - Tr,[exp{-f]H,,,,)]} . (6) 

Tr2 denotes the trace over all two neutron states, H is the full Hamiltonian, and -fffree is 
the free Hamiltonian. By integrating over the center of mass momentum and enforcing 
spherical boundary conditions on the relative displacement between the two particles, the 



density of scattering states can be related to the total elastic phase shift S{E) [35|, 



62(T)-6f^(T) 



/9 



2V2 



dE e ^^^'^6{E) + bound state contribution. 



(7) 



If there are two-body bound states in the spectrum with binding energies Es^i, there is an 
additional contribution 



2V 



1/2 



,/3 Ss, 



1 



(8) 



In the unitary limit, where the effective range is zero and scattering length is infinite, we get 



62(T) = 3x2"^^ 0.530. 



(9) 



For zero effective range but arbitrary scattering length Cscatt the second virial coefficient is 

fj[l-erf(|x|)]-^ forx<0, 

^[l-erf(x)] forx>0, 

where erf is the error function, Eb is the two-particle binding energy for positive scattering 
length, and 

x= . (11) 



(12) 



2vrascatt 

As the effective range goes to zero we have the relation 



\E 



1 



B\ 



ma. 



2 ' 

scatt 



and therefore we can write 



5j[l-erf(|x|)]-^ forx<0 



4^2 



^[l + erf(x)] - ^ forx>0. 



(13) 



III. ONE-DIMENSIONAL MODEL 



We begin our analysis of finite cutoff errors with a one-dimensional model. The model 
consists of nonrelativistic spin-1/2 fermions with an attractive zero-range interaction and is 
the continuum limit of the attractive one-dimensional Hubbard model. Both the attractive 
and repulsive versions of the one-dimensional Hubbard model have been studied in the 



literature 
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42l |. We consider the attractive case as a toy model 



of short-range attractive forces in nuclear systems. As we will see, the problem of large 
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FIG. 1: The connected amputated two-particle Green's function. 

discretization errors appears even in this one- dimensional model which has no ultraviolet 
divergences. 

In the continuum limit the Hamiltonian has the form 

H = — / dx al{x)-—^ai{x) + C / rfa; a|(a;)a|(a;)a-f(a;)aj(a;), (14) 

where m is the mass, C < 0, and and a| are annihilation and creation operators for spin i. 
The connected amputated two-particle Green's function equals the sum of bubble diagrams 
shown in Fig. [H Any connected scattering process consists of two-particle Green's functions 
linked together with free particle propagators. 

Let G2{po,Px) be the amplitude for the connected amputated two-particle Green's func- 
tion, where po is the total energy and px is the total spatial momentum of the two particles. 
We sum the bubble diagrams in Fig. [T] and get 

—iC 1 
1-tC - Uipo^Px) + U{po,Px) 



where 



, dqodqx i i . . 

n(po,Px) = / -—T 7— -^2 X — . (16) 



In the continuum limit we find 



n(po,Px) = — , ^ , (17) 

G2(P0,P.) = 1 ^ m ■ (18) 
Since C < there is a bound state pole in the two-particle Green's function at energy 
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We can obtain the same result by solving the Schrodinger equation for the two-particle 
system in the center of mass frame. If x is the relative separation between particles then, 
in the center of mass frame, the reduced Hamiltonian is 

1 



HcM 

and the ground state wavefunction is 



m dx^ 



+ C5{x), 



oc exp ( -^C* \x\ 



The ground state energy with center of mass kinetic energy included is 



Eo{Px) = 



4 4m' 



(20) 



(21) 



(22) 



where Px is the total momentum. 

The one-dimensional system is finite in the continuum limit and therefore no regulariza- 
tion nor renormalization is needed. Nevertheless we impose an ultraviolet cutoff on the 
momentum in order to study the resulting cutoff errors. With a momentum cutoff at A we 
find 



n(po,Pa;,A) 



X 



(23) 



where is some homogeneous combination of the parameters mpo and p^. The combination 
will depend on the details of the chosen regularization scheme. The regularized two-particle 
Green's function has the form 

1 



G2(P0,Px, A) 



iC{A) 



X 



1 + % 



(24) 



2i/mpo-^ 

We define the scale-dependent couphng C(A) so that the pole in the rest frame remains 
exactly at 

Po = ^- (25) 



A. One- and two-particle lattice dispersion relation in one dimension 

We investigate the cutoff error in more detail using a Hamiltonian lattice formalism. On 
the lattice the cutoff momentum scale A corresponds with Tra~^, where a is the lattice spac- 
ing. Throughout our discussion of the lattice formalism we use dimensionless parameters 
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and operators, which correspond with physical values multiplied by the appropriate power 
of a. However final results are reported into physical units. We start with the simplest 
possible lattice Hamiltonian giving (fT^ in the continuum limit. We let 

i7(o) = i^(o) + 0o)^ (26) 
K^^^ = —y^a\{n^)ai{n^) - -^y^ \a\{n^)ai{n^ , (27) 

^(0) = c5^a|(nJa{(nJatK)aiK). (28) 

We refer to H^^^ as the standard lattice Hamiltonian. The zero superscript signifies that 
it is the simplest possible lattice formulation. Later in our discussion we consider more 
complicated lattice actions. We choose the mass to be m = 939 MeV and fix the lattice 
spacing at a = (50 MeV)^^. This corresponds with a cutoff momentum of A = na^^ 157 
MeV. 

The single-particle dispersion relation for the standard lattice Hamiltonian is given by 

J^)(p^) = 1 X (1 - cos p.) = ^ + 0{pt). (29) 
m 2m 

In Fig. [2]we have plotted ij^^\px) versus the continuum result uj{px) = Px/ (2m) for momenta 
in the first Brillouin zone \px\ < A. The relative error between cu^^^ and cu is 10% or less for 
\Px\ < A/3. 
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FIG. 3: The lowest two-particle energies for the standard lattice action, and Ei , and the 
corresponding continuum limit values, Eq and Ei. The continuum coupling is C = —0.0400 while 
the lattice couphng is C(A) = -0.0407. 



We now consider two-particle states with one spin-up particle and one spin-down particle. 
We start with a small value for the coupling, C = —0.0400. Let be the total momentum 
of the two-particle system. In the continuum limit the ground state energy at zero total 
momentum is 

Eo{p^ = 0) = — = -0.3756 MeV. (30) 

It is convenient to place the two-particle system in a periodic box. We choose the box 
length to be L = 1 MeV~^. Since the ground state wavefunction depends on the relative 
separation x as 

ipoix) (X exp (^"iT^C \x\^ ^ exp [—(19 MeV) ■ |x|] , (31) 

the effect of the boundary at L = 1 MeV"^ on the ground state energy is negligible. The 
box length does however determine the level spacing between unbound scattering states. 

For the lattice calculation the interaction strength is tuned so that the ground state 
energy in the rest frame matches the continuum result of —0.3756 MeV. This gives an 
adjusted coefficient of C(A) ~ —0.0407. In Fig. [3] we show the two lowest energy levels 
of the two-particle system as functions of px- The two-particle ground state and lowest 
scattering state energies for the lattice are labelled E^^^ and Ef'^ respectively, while the 
corresponding continuum limit values are labelled Eq and Ei. The mismatch between lattice 
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FIG. 4: The lowest two-particle energies for the standard lattice action, and Ei , and the 
corresponding continuum limit values, Eq and Ei. The continuum coupling is C = —0.1000 while 
the lattice coupling is C(A) = —0.1105. 



and continuum results for the two-particle energies is roughly the same size as the mismatch 
between single-particle lattice and continuum kinetic energies, cu^'^^ and cu. 

Keeping other parameters the same we now repeat the two-particle energy calculations 
at stronger coupling, C = —0.1000. In this case the continuum limit ground state energy 
at Px = is 

Eo{p^ = 0) = — = -2.348 MeV. (32) 

Tuning the lattice interaction to match this ground state energy gives an adjusted coefficient 
of C(A) ^ —0.1105. Results for the two-particle ground state and lowest scattering state 
are shown in Fig. |H While the agreement for the ffist excited states e[^^ and Ei has not 
changed noticeably, the deviation between lattice and continuum results for the ground state 
energy has increased substantially for px 7^ 0. We have chosen the lattice coupling C(A) so 
that E^^ = Eq at p^ = 0, and so the disagreement between E^^ and Eq is proportional p^. 
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B. Broken Galilean invar iance 



In the continuum limit Galilean invariance requires that the total energy of the two- 
particle system rises quadratically with the total momentum p^, 

E{p^) = E{0) + ^. (33) 

However any regularization scheme with a preferred reference frame breaks Galilean invari- 
ance to some extent. In the following we show how broken Galilean invariance on the lattice 
can result in large cutoff errors at strong coupling. 

Let the momenta of the two particles be Px/'i — Qx and Px/'2 + Qx, where Px is the total 
momentum and 2qx is the relative momentum between the particles. We consider first the 
ground state. The average value of the relative momentum 2qx for the two-body ground 
state wavefunction grows proportionally with m \C\. For C = —0.0400 we have m \C\ = 37.6 
MeV, while for C = -0.1000 we have m \C\ = 93.9 MeV. For the latter case is not 
small compared with the cutoff momentum 157 MeV. It is not so large as to invalidate the 
assumption that we have a sensible low-energy effective field theory. However it is large 
enough that one of the constituent particle momenta can reach the Brillouin zone boundary 
at ±A even though \px/2\ is less than A. At the zone boundary the lattice kinetic energy 
tu^^'^ is significantly lower than the continuum kinetic energy tu. This error in the dispersion 
relation produces a ground state energy E^\px) which is lower than the continuum result 
Eq[Px) at strong coupling. This is the effect we observe in Fig. HI 

The problem with large relative momentum does not occur for low-energy scattering 
states above the ground state. This is because the wavefunctions for these scattering states 
are peaked around the asymptotic momenta of the particles, Px/2 — Qx and Px/'i + Qx- In 
the infinite L limit we have 

P 1 fP^ 1 fPx , Y pI ,ql (^.s 

ql = mE- ^. (35) 

If mE and p^. are much less than A^, it follows that is also much less than A^. Hence 
the single-particle momenta Px/2 — qx and Px/2 + qx are small compared with A, and cutoff 
errors should remain small for low-energy scattering states even at strong coupling. 
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FIG. 5: The two-particle energies for the lattice, e'^^ and ^ and continuum, Eq and E\. The 
continuum coupling is C = —0.1000 while the lattice coupling C(A) is set by matching Eq{px) at 
\p^\ =p™« ~ 60 MeV. 



C. Cutoff errors at nonzero temperature 

In the classical regime the equipartition theorem tells us that the distribution of momenta 
satisfies 



n2 



1, 



= pT = (36) 

Since the cutoff errors are proportional to p^, this suggests that at fixed lattice spacing 
the cutoff errors for the dilute system should increase linearly with T. One approach to 
removing this error at nonzero T is to define the lattice coupling C(A) by matching the 



continuum energy Eq{px) at \px\ = \JlmT rather than at px = 0. In Fig. [5] we show the 
two-particle energies when C(A) is fit to Eq{px) at \px\ = p™^ — 60 MeV. This redefinition 
has the unattractive feature that the lattice coupling C(A) is now temperature dependent. 
Furthermore it does not fix the problem of strongly-broken Galilean invariance. The cutoff 
error has simply been shifted to momenta \px\ 7^ P™- However it does remove large cutoff 
errors from the lattice simulation at temperature T. This is essentially the approach used 
in 0, Is], where the lattice coupling C(A) was determined by matching the continuum limit 
value for the second virial coefficient b2{T), where T is the chosen simulation temperature. 

There are other techniques which actually reduce the breaking of Galilean invariance 
on the lattice. One possibility is to remove all of the O j dependence using higher- 
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dimensional operators. This includes two- derivative interactions such as 



d 



aj(x)aj(x) — [a|(x)a|(x)] 



(37) 



O^ij) = a|(x)a|(x) 
(i^a|(x) 



da^{x) da^^x) 



dx dx 



+ at(x 



dx'^ 



+ 



dx"^ 



a' X - 2 



da}^{x) (ia|(x) 



. d'^a\[x) 

+ a (x)— V 
dx'^ 



a^{x)a^{x). 



(38) 



Oi could be tuned to cancel the broken Galilean invariance while O2 could be tuned to 
reset the effective range to zero. However new interactions such as these can introduce sign 
oscillations and other complications in Monte Carlo simulations. Therefore we first try a 
less expensive approach where the interaction is left alone and only the lattice kinetic energy 
is modified. 



D. 0(a^ )-improved and 0(a^)- well-tempered actions in one dimension 

Let us consider replacing the standard lattice kinetic energy action in (1271) with an O(a^) 
improved kinetic energy with next-to-nearest neighbor hopping, 

1 



^ X — ^^a\{n^)ai{n^) - | x 



4 



2m 



\ a\{n^)ai{n:, + 1) + a\{n^)ai{n^ - 1) 



I2 ^ 2m ^ [al(^2:)«i("'x + 2) + a\{nx)ai{n^ - 2) 



This gives the dispersion relation 

1 



(39) 



m 



5 4 1 ,0 

- - -cosp:, + 12 "^^^ 



2m 



+ o {pD . 



(40) 



Matching Eq{px = 0) = —2.348 MeV for the improved lattice action gives an adjusted coef- 
ficient of C(A) = —0.1173. Results for the two-particle ground state and lowest scattering 
state with the improved action are shown in Fig. [61 We see that the deviation between 
lattice and continuum results for the ground state energy has been reduced for Px ^ 0- 

While the errors for the improved action are better than that for the standard action, 
better agreement seems possible. Instead of removing the O (p^.) term from the lattice 
dispersion relation, this time we tune the coefficient of the O (p^) term to match as best as 
possible the continuum dispersion relation uj{px) = pi/ (2m) over the full range —A < Px ^ ^■ 
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FIG. 6: The lowest two-particle energies for the improved lattice action, Eq and E\ , and the 
continuum energies, Eq and Ei. The continuum coupling is C = —0.1000 while the lattice coupling 
is C(A) = -0.1173. 



Let us define the 0(a^)-well-tempered lattice kinetic energy action K^""^^^ and dispersion 
relation uj^'"^'^^ (p^) , 

j^i^ti) ^ ^(0) ^ ^ (^(1) _ ^(0)) ^ (41) 

= + s {u;^'\p,) - , (42) 

where the unknown coefficient s is given by the integral constraint 

A A 

jdp,u^-''\p^) = jdp.^. (43) 

-A -A 

Solving for s gives s = |7r^ — 4 ^ 2.5797. We show a comparison of the lattice dispersion 
relations uj^^\ uj^^\ u;^™*^-', and the continuum limit cu in Fig. [71 Matching Eo(p^ = 0) = 
—2.348 MeV for the well-tempered lattice action gives an adjusted coefficient of C(A) = 
—0.1260. Results for the two-particle ground state and lowest scattering state for the well- 
tempered action are shown in Fig. [HI The deviation between lattice and continuum results 
for the ground state energy has been substantially reduced. 
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FIG. 7: Comparsion of the lattice dispersion relations uj^^\ a;^^*-^^ and the continuum limit 
(jj. 
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FIG. 8: The lowest two-particle energies for the well-tempered lattice action, ' and E\ ' , 
and continuum, Eq and Ei. The continuum coupling is C = —0.1000 while the lattice coupling is 
C(A) = -0.1260. 

IV. ZERO-RANGE NEUTRONS IN THREE DIMENSIONS 

We now explore how various lattice actions affect cutoff errors in three dimensions. We 
consider spin-1/2 fermions in three dimensions with zero- range attraction. This gives an 
approximate description of interacting neutrons below 1% of nuclear matter density. To 
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demonstrate the generality of our lattice error analysis we consider both the grand canon- 
ical ensemble in the Euclidean lattice formalism as well as the canonical ensemble in the 
Hamiltonian lattice formalism. 

In the continuum limit the Hamiltonian for zero-range neutrons has the form 

iJ = --L^ f d^ral{f)V^ai{r) + C j d^r a\{7^a\{r)a^{f)ai{r). (44) 

Just as in our one-dimensional model, the connected amputated two-particle Green's func- 
tion for zero-range neutrons in three dimensions is given by the sum of chained bubble 
diagrams shown in Fig. [TJ Any connected scattering process can be constructed from 
two-particle Green's functions linked together with free particle propagators. While the 
two-particle Green's function is divergent, all of the new loop integrations produced by 
connecting two-neutron Green's functions are ultraviolet finite. 

Let G2{po,p) be the amplitude for the connected amputated two-particle Green's function, 
where po is the total energy and p is the total spatial momentum of the two particles. We 
sum the bubble diagrams in Fig. [T] and get 

—iC 1 

G2{Po,P) = z TTTTT ^ = 1 , rtr ^' (^^) 

where 

n (Po.Pl = / X i-^^. (46) 

Since Il{po,p) is ultraviolet divergent we renormalize the coupling C to absorb the diver- 
gence. In the end we get 

G2iPo,P,A) = ^^i^ — (47) 



where agcatt is the s-wave scattering length and is some homogeneous combination of the 
parameters mpo and p^ which depends on the regularization scheme. The cutoff error can be 
regarded as a momentum/energy-dependent 0{Q'^/A) modification to the inverse scattering 
length. 
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A. One- and two-particle lattice dispersion relation in three dimensions 



Just as in the one- dimensional case we start with the simplest possible lattice Hamiltonian 
that reproduces (jUj) in the continuum limit. We let 

i/(o) = K(o) + v^(o)^ (48) 
= ~ X] aj(ns)aj(n^) - ^ ^a\{ns)ai{ns + Is) + a\{ns)ai{ns - Is) , (49) 
y(o) =cY,(^\ins)a\{ns)a^{ns)ai{ns). (50) 

Here is a three-dimensional spatial lattice vector and Is = x,y,z are lattice unit vectors 
in each of the 3 spatial directions. The s subscript is our notation for spatial lattice vectors 
with no time component. H^^^ is the standard lattice Hamiltonian. The zero superscript 
again signifies that it is the simplest possible lattice formulation. We take the same values 
m = 939 MeV for the neutron mass and a = (50 MeV)~^ for the lattice spacing. This 
again yields A = ira^^ ^ 157 MeV for the cutoff momentum. The single-particle dispersion 
relation for the standard lattice Hamiltonian is 

^'''^P^^ = i E (l-eospJ = |^ + 0(b-n. (51) 

l3=x,y,z 

We consider a lattice system which is a periodic cubic lattice of length L. If we set 
the two-body scattering pole in the rest frame at energy -Epoie then the cutoff-dependent 
coefficient C(A) satisfies 

1 1 1 



lim 



C{A) L^oo L3 ^ _Epole + 2^(0) i27Tks/L) ' 

ks 

where the components of ks are integers from 0,1,---,L — 1. If there is a two-body bound 
state then we can take Ep^ie equal to negative the binding energy. Alternatively we can 
choose -Epoie to be the pole nearest threshold and use Liischer's formula for the finite volume 



two-body spectrum 



43 



44| 



7-1 / rN 47rascatt r-i C^scatt '^scatt ,1 frn\ 

ML) = —-^[1 - c,-^ + c,-^ + ■■■], (53) 



where Ci = -2.837297, C2 = 6.375183. 
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FIG. 9: The lowest two-particle energies for the standard lattice action, E^^ and and the 

corresponding continuum limit values, Eq and Ei. 

As an example we set E^^ic = -0.300 MeV and find that C(A) = -9.464 x 10"^ MeV'^. 
This corresponds with a scattering length flgcatt = 11-76 fm. For L = 30 we plot the energy 
of the lowest two energy states E^^ and Ef^ as a function of the total momentum in Fig. [9] 
and compare with the corresponding continuum limit values, Eq and Ei. In physical units 
L = 30 corresponds with 118 fm. This is about ten times the scattering length and so finite 
volume effects are negligible. The box length does however determine the level spacing 
between unbound scattering states. We take the total momentum pg along the x-axis so 
that ps = {px, 0, 0). Just as we found in the one- dimensional model at strong coupling, we 
encounter the same problem of broken Galilean invariance. While the agreement between 
the excited states Ef''^ and Ei is not bad, the deviation between lattice and continuum 
results for the ground state energy is substantial for p^ ^ 0. Since we have chosen the 
lattice coupling C(A) so that Eq^^ = Eq at p^ = 0, the disagreement between E^^ and Eq is 
proportional p^.. 
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B. 0(a^ )-improved and 0(a^)- well-tempered actions in three dimensions 



Just as in the one- dimensional case we can replace the standard lattice kinetic energy 
with an 0(a^)-improved kinetic energy, 



aUns)ai(ns) - 



3 2m ^ I ^ 



Am 

+ ^ X ^ ^ ^a\{ns)ai{ns + 2/^) + a\{ns)ai{ns - 2/^) 

This gives the dispersion relation 

,(1) 



{fis + Is) + a\{ns)ai{ns - Is) 



(54) 



5 4 1 , , 

- - - cos + Y^cos (2piJ 



\ps) = - E 
m ^-^ 

ls=x,y,z 

In this case the renormalization condition for C(A) is 



|;+o(ip-r). 



(55) 



= lim — \ 

C(A) L^oo L3 ^ _E^^^^ + 2ujW{2^kjL) ' 



(56) 



and for Epoic = -0.300 MeV we find C(A) = -1.1031 x 10"^ MeV"^. Resuhs for the 
two-particle ground state and lowest scattering state with the improved action are shown in 
Fig. [ini The results are somewhat better for the 0(a^)-improved kinetic energy, though the 
agreement between E^^ and Ei all the way up to the cutoff momentum should be regarded 



19 




as accidental. As in the one-dimensional case, we expect that better agreement may be 
possible for the ground state using a well-tempered action. 
We define the 0(a^)-well-tempered kinetic energy as 

where s is given by the following integral constraint on the resulting dispersion relation: 

AAA AAA ^ 

j j j dp^dpydp, ^-"''^1),) = j j jdp^dpydp,^. (58) 

-A-A-A -A-A-A 

Since both uj^'^^^\ps) and p^/{2m) decompose as a sum of separate terms for p^, py, and Pz, 
this gives the same result as in the one-dimensional case, s = |7r^ — 4 ^ 2.5797. For the 
well-tempered action 



lim 



1 >^ 1 

C(A) £Z^;.L^ -Epoie + 2a;(-")(27rfc,/L)' ^ 

ka integer 

and for Ep^ie = -0.300 MeV we find C(A) = -1.3273 x 10"^ MeV'^. Results for the two- 
particle ground state and lowest scattering state with the well-tempered action are shown 
in Fig. [m Just as in the one-dimensional model, we find the deviation between lattice and 
continuum results for the ground state energy has been substantially reduced. 
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The well-tempered kinetic energy appears to fix much of the strongly-broken Galilean 
invariance on the lattice. In the remainder of our analysis we see if it also fixes the large 
discretization errors at nonzero temperature. We do this by calculating the second virial 
coefficient 62(^)5 which was found to have large discretization errors in ^. There are 
several ways to calculate b2{T) on the lattice, and it is not obvious that the lattice errors 
are the same for different calculations. Therefore in the next two sections we consider two 
different lattice calculations of b2{T). The first calculation relies on the virial expansion of 
the density in the grand canonical ensemble, 

p=^[z + 2h{T)z^ + ---]. (60) 

We use the Euclidean lattice formalism with nonzero temporal lattice spacing for this 
calculation. The second method finds b2{T) by means of the two-particle partition function, 

b2{T) - b'riT) = ^ {Tr2[exp{-[3H)] - Tr2[exp{-f3H,,,,)]} . (61) 

We use the Hamiltonian lattice formalism for this calculation. 



V. GRAND CANONICAL EUCLIDEAN LATTICE CALCULATION FOR 62 (T) 

In this section we calculate the second virial coefficient b2{T) in the grand canonical 
ensemble using the Euclidean lattice formalism. We use the same values m = 939 MeV for 
the neutron mass and a = (50 MeV)~^ for the lattice spacing. We set the temporal lattice 

n n 

spacing at at = (24 MeV) ^. These are the same values as used in [7|, l8|. We define at 
as the ratio of temporal to spatial lattice spacings. In our notation n = {rit, Us) denotes 
space-time lattice vectors, c and c* are Grassmann variables for the neutrons in the path 
integral formalism. is a lattice unit vector in the temporal direction. Ig = x,y,z are 
lattice unit vectors for the spatial directions. Also fi is the chemical potential, L is the 
spatial length of the cubic lattice, and Lt is the temporal length. 

In the grand canonical ensemble the partition function can be written as 

Z (X J DcDc* exp (-5) , (62) 

where 

S = 5frec + C'ate^^"^^ c\{n)c\{n)c^{n)c^{n), (63) 

n 
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and the standard free lattice action is given by 



sSL = E bnn)c^{fi + 6) - e'^-c*(n)Q(n)] + a,e>^^^ x - 



m 

n.i 



L s [<=•■' 



n]Ci{n + Is) + c*(n)ci(n - Is 



n,la,i 



Our coupling constant C differs from the coupling constant C appearing in (4, |7|, l8| 
ever the two are simply related, 

2 



Let us define 



Then we have 



5'static = [c*(n)ci(n + 0) - e'^"'c*(n)ci(n)] . 



c(0) _ c . _|_ ^ pM"t c(0) 
'-'free ~ static + ^kinetic 



where 



0(0) 
kinetic 



n cAn - 



2m 



c*(ra)Q(ra + /s) + c*{n)ci{n - Is) 



'5'ki'netic is analogous with i^*^*^) in (H^ . The 0(a^)-improved action has the form 



"^free " static + 0*6 "JkinetiC 



where 



9^') - - X 

"-^kinetic ■ ^ 



! - ^ yic*(n)Q(n) - ^ X 
4 m 3 2m 



c*(n)ci(n + ls) + c*(n)ci(n - Is 



+ ^ X ^ 5^ |c*(^)ci(^ + 2/.) + c*(n)Q(n-2/,) 



The 0(a^)-well-tempered action has the form 



Q(wtl) _ Q I ^, ^iiat c(wtl) 

'-'free " 'Jgtatic "t" "tc "^kinetic 



where 



c.(wtl) ^ c.(0) / _ c^(0) ^ 

kinetic kinetic ' ^Vivir-fir- 



kinetic kinetic / ' 
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s = -vr^ - 4 ^ 2.5797. 



(72) 



For each lattice action we define the free neutron propagator, 

n J DcDc* di{k)c*{-k)exp{-S{,^,) 

Dfreeik) = fn n * / c ^ ('^o sum on t), 

J DcDc* exp (-5fi.cc) 



(73) 



where the components of A; = (/eq, kg) are integers. Our conventions for the lattice Fourier 
transform are 



^•^"tfep ■27rns-fcs 



~f{k) = Y,e^e 



'fin) 



fin) 



LtL- 



e 



-fik). 



(74) 
(75) 



Let uj{27rks/ L) be the lattice dispersion relation, either cu^^\2Trks/ L), ij^^\2nks/ L), or 
Lj(wti)(27rA;s/-L) as defined in the previous section. Then we have 

1 



Direeik) 



a 



e * Lt _ (.fiat atet"^^uj{27rks/L) 



(76) 



In it was shown that the cutoff-dependent coupling constant is given by the constraint 



atC'{A) L^oc L3 



lim -^y^ 



-at-Epolo 



1 — atU!{2nks/L) 



1 2 • 



(77) 



We use the Euclidean lattice action to compute the neutron density as a function of 
temperature, chemical potential, and interaction strength. Let p^ee be the free neutron 
density and p be the neutron density with interactions. Then from the virial expansion we 

get 

P - Pfree = ^ HT) - ^^''(T)] + 0{z'). (78) 

We note that our convention for the second lattice coefficient &2(r) is slightly different from 
the one used in [tI. The densities pfree and p are computed using the free and full neutron 
propagators. 



1 d 



LtL^ 



^ Dirceik)e 



l3L^dp 



\nZ = 2 



e ^« 



(79) 
(80) 
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n = 



1 2 



n 



FIG. 12: Two-particle bubble diagrams contributing to the neutron self-energy to order z'^ 



The full neutron propagator D{k) can be expressed in terms of the neutron self-energy, S(A;), 



D{k) 



(81) 



1 - S(A;)Dfree(fc) 

We compute the self-energy to order by summing the two-particle bubble diagrams shown 
in Fig. [121 Further details of the calculation can be found in \\. 

In addition to these local actions we also consider a dispersion relation given by 

(82) 



ls=x,y,z 



where 



kl=mod{ki^,L), Ik'A <L/2. 



(83) 



This quadratic dispersion relation was used in jlSj to reduce cutoff effects. Since it equals the 
continuum dispersion relation for 27[ks/L < A we expect it also to be effective in reducing 
errors due to broken Galilean invariance. However it does not correspond with a local lattice 
action. It was implemented in [l3| by Fourier transforming back and forth between position 
space and momentum space. Unfortunately this results in a steeper computational scaling 
for the Monte Carlo algorithm as a function of volume. Nevertheless there is no significant 
computational problem for the perturbative calculation presented here, and so we include it 
in our analysis for comparison. 

We can compute b2(T) at any small fugacity, and so we choose z = ~ 0.0067. We take 
the lattice length to be L = 8, which is sufficiently large enough that the finite volume error 
for the local actions is less than 1%. The non-local action associated with cu^'^"^'^) appears 
to have a slightly larger finite volume error. Using each of these dispersion relations, we 
compute the second virial coefficient for a range of scattering lengths, Ogcatt = —4.675, —9.35, 
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FIG. 13: Plot of b2iT) for uj^°\ w^^"^'^) and the continuum limit for T = 1.0 MeV. 
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FIG. 14: Plot of 62(r) for cj(™*^), ^^(q"^^) ^.^^ continuum limit for T = 2.0 MeV. 



— 18.70, +18.70, +9.35 fm. In Fig. [T3]we show the resuhs for b2{T) as a function of inverse 
scattering length for dispersion relations uj^^\ uj^^\ iu^'^^^\ and uj^'^^^'^^ at T = 1.0 MeV. We 
also show the continuum limit result given in (fT3l) . Analogous results at temperature T = 
2.0 MeV are shown in Fig. [HI We see that of the various lattice dispersion relations, tu*^™*^^ 
and uj^'^^^'^^ come closest to the continuum limit. 
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0.5 




r(MeV) 

FIG. 15: Plot of Aog^^^-j as a function of temperature for the dispersion relations u!^^\ ui^^\ 
and w^-^"^^). 



We can compare the different lattice actions in a slightly different way. Let us think of the 
temperature as fixed and the scattering length as varying. When deriving (H7ll we found 
that the finite cutoff error can be regarded as a momentum/energy-dependent 0{Q'^/A) 
modification to the inverse scattering length. In the continuum limit b2iT) = 3 x 2^2 at 
infinite scattering length for all T. At finite cutoff let a^att(^) tie the scattering length 
for which 62(T) = 3 x 2-2 ^ 0.530. We can now interpret the cutoff error as a small 
modification to the inverse scattering length, Aa~att(^) = "l/^^^att- ^^e continuum 
limit Aascatt(^) = T, and the shift provides a simple quantitative measure of the 

cutoff error near infinite scattering length. 

We expect broken Galilean invariance due to the cutoff to introduce a term of size 0(pf /A) 
in Aa~^]^^^{T). In the classical regime we know from the equipartition theorem that the 
average value of scales linearly with the temperature T. Therefore we expect Aa'^^tl^) 
also to scale linearly with T. In Fig. [15]we plot Aa~f.\^^{T) for the dispersion relations u^^\ 

see the expected linear dependence in Aag(,att(^) for small T. 
We also see that Aa-^i^^{T) for cu^™") and are quite a bit smaller that Aa^^^^^iT) for 

u''^^ and u;'-^^ In fact most of the cutoff error at nonzero T appears to have been removed. 
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FIG. 16: Plot of 62 at T = 1.0 MeV for the standard Hamiltonian lattice action at T = 1.0 MeV. 
For comparison we show results for the standard Euclidean lattice action, the standard Hamiltonian 
lattice action with Galilean invariance imposed by hand, and the continuum limit. 



VI. TWO-PARTICLE HAMILTONIAN LATTICE CALCULATION FOR 62 (T) 

In this section we return to the Hamiltonian lattice formalism and compute fe2(7') using 
the two-particle trace formula, 

62(T) - 6f-(T) = ^ {Tr2[exp(-/3i7)] - Tr2[exp(-/3i/f,ee)]} • (84) 

As before we take a = (50 MeV)~^ and L = 8. In Fig. [16] we show results for the 
standard Hamiltonian lattice action at T = 1.0 MeV. For comparison we show results for 
the standard Euclidean lattice action, the continuum limit, and the standard Hamiltonian 
lattice action with Galilean invariance imposed by hand. We impose Galilean invariance 
by computing the spectrum of H in the rest frame. We then boost the result for nonzero 
total momentum Ps using 

E{p,) = E{0) + ^. (85) 

We see that both the standard Hamiltonian lattice results and standard Euclidean lattice 
results deviate from the continuum limit by about the same amount. We also see that the 
standard Hamiltonian lattice action with Galilean invariance is almost identical with the 
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FIG. 17: Plot of 62 at T = 1.0 MeV for the well-tempered Hamiltonian lattice action at T = 1.0 
MeV. For comparison we show results for the well-tempered Euclidean lattice action, the well- 
tempered Hamiltonian lattice action with Galilean invariance imposed by hand, and the continuum 
limit. 

continuum limit. This suggests that broken Gahlean invariance is in fact responsible for 
most of the cutoff error at T = 1.0 MeV. 

We show in Fig. [Tfl results for the O (a^)-well-tempered action at T = 1.0 MeV. The 
four curves shown are for the O (a^)-well-tempered Hamiltonian action, O (a^)-well-tempered 
Euclidean lattice action, the continuum limit, and the O (a^)-well-tempered Hamiltonian 
lattice action with Galilean invariance imposed by hand. In this case all four curves agree 
rather well. The well-tempered action clearly preserves Galilean invariance much better 
than the standard action and removes most of the cutoff error in both the Hamiltonian and 
Euclidean lattice formalisms. 

VII. SUMMARY AND DISCUSSION 

In this study we investigated the temperature dependence of lattice discretization errors 
in nuclear lattice simulations. As a warm-up exercise we started with the one-dimensional 
attractive Hubbard model. We found that when the interaction was strongly attractive the 
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dispersion relation for the two-particle ground state showed significant cutoff errors. This 
cutoff error could be attributed to the breaking of Galilean invariance. The same problem of 
strongly-broken Galilean invariance was found in three dimensions for interacting neutrons 
with an attractive zero-range potential. 

We showed that part of the error due to broken Galilean invariance could be eliminated 
by using an 0(a^)-improved lattice kinetic energy. The 0(a^)-improved action includes 
next-to-nearest neighbor hopping terms in order to match the single particle dispersion 
relation pi/ (2m) up to terms O While the improved action was better than the 

standard action, we found even better results when using an 0(a^ )-well-tempered kinetic 
energy lattice action. The 0(a^)-well-tempered action includes the same next-to- nearest 
neighbor hopping terms as the 0(a^)-improved action. However in this case the coefficients 
of the various terms are adjusted to match the integral of pf/(2m) over all momenta below 
the cutoff, 

AAA AAA ^ 

j j j dp^dpydp, J^''\ps) = J J Jdp^dpydp,^. (86) 

-A-A-A -A-A-A 

We then performed two separate calculations of the second virial coefficient b2{T) using 
the various different lattice actions. In the first calculation we extracted b2{T) in the grand 
canonical Euclidean lattice formalism using the virial expansion of the density. In the second 
calculation we determined b2{T) by a Hamiltonian lattice calculation of the two-particle 
partition function. In both cases we found that the 0(a^)-well-tempered lattice action was 
superior to both the standard action and 0(a^)-improved lattice action. In fact the well- 
tempered action reduced the temperature-dependent cutoff errors as much as the non-local 
action favored in [isl. This non-local action corresponds with the quadratic dispersion 
relation uj^'^^^'^\ However the 0(a^)-well-tempered lattice action has the advantage of being 
a local action. Therefore it can be implemented in most Monte Carlo lattice algorithms 
without increasing the computational scaling. 

The well-tempered action is a simple way to reduce lattice errors at nonzero temperature. 
While the discussion here has focused on zero-range pionless effective field theory, it seems 
clear that the well-tempered action fixes the problem of strongly-broken Galilean invariance 
quite generally. With this increase in accuracy it should now be possible to perform an 
accurate lattice calculation of the third virial coefficient b^lT). b?j{T) was recently calculated 



for two-component fermions in limit of zero effective range and infinite scattering length 
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and this calculation can now be checked on the lattice. 

We can organize the various kinetic energy actions introduced here somewhat more sys- 
tematically. In the Hamiltonian lattice formalism let 

Kj-hop XI ^al{ns)ai{ns + jig) + a\{ns)ai{ns - fh) (87) 

for integers j > 0. In the same way in the Euclidean lattice formalism let 
Sj-hop = ^ XI ['^i (^)'^i(^ + ^L) + c*{n)ci{n - jis) 



(88) 



for integers j > 0. For any chosen lattice action K^^^ or ^S^^ we assign a set of hopping 



coefficients vj"^'* such that 



or 



i=0,l,2,- 



'S'kinetic ~ X ^ ^y'"\ ^ Sj. hop- 
i=0,l,2,- 

Then the corresponding single-particle dispersion relation is 

1 



7 ("■) 
■'Vj ' COS 



(89) 



(90) 



(91) 



jr=0,l,2,--- ls=x,y,z 

The 0(a^)-improved action is defined so that its dispersion relation u^'^\ps) agrees with 
01/ {2m) up to terms 0[\Psf). The 0(a^)-well-tempered action is defined so that its 
dispersion relation uj'''"^'^\ps) satisfies 

u^'^'^\2T:k/L)=J''\2T:k/L)+u^^\2T:k/L) + s\u^^\2T:k/L)-J^\2T:k/L^ 

AAA AAA _^ 

j j jdp,dpydp,UJ^^''\ps)^ j j jdp^dpydp,^. (93) 

-A-A-A -A-A-A 

The generahzation to higher-order actions is straightforward. The hopping coefficients for 
the various actions up to 0{a^) are shown in Table 1. 
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standard 


0(a^)-improved 


0(a^)-well-tempered 


0(0"^) -improved 


0(a^)-well-tempered 




1 


5 
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6 


49 

36 
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Vl 
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4 
3 


27r2 1 
9 3 


3 
2 


TT^ 13 
4 24 


V2 





1 

12 


7r2 1 
18 3 


3 

20 


TT^ 2 

10 3 


V3 











1 

90 


7r2 1 
60 8 



Table 1: Hopping coefficients for kinetic energy lattice actions up to 0{a^) 
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We note that vq for the well-tempered action is 7r^/6 for all orders. This is because the 
integral of cos {2t: jki J L) vanishes for j 7^ and so only the vq term survives in the integral 
over momenta. These higher-order well-tempered actions may be useful if we wish to use the 

same lattice action for high-accuracy nucleon-nucleon scattering phase shifts and many-body 
simulations at nonzero temperature. 
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